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ABSTRACT 

The power density spectrum of a light curve is often calculated as the average of 
a number of spectra derived on individual time intervals the light curve is divided 
into. This procedure implicitly assumes that each time interval is a different sample 
function of the same stochastic ergodic process. While this assumption can be applied 
to many astrophysical sources, there remains a class of transient, highly nonstationary 
and short-lived events, such as gamma-ray bursts, for which this approach is often 
inadequate. The power spectrum statistics of a constant signal affected by statistical 
(Poisson) noise is known to be a in the Leahy normalisation. However, this is no 
more the case when a nonstationary signal is also present. As a consequence, the 
uncertainties on the power spectrum cannot be calculated based on the properties, 
as assumed by tools such as XRONOS powspec. We generalise the result in the case 
of a nonstationary signal affected by uncorrelated white noise and show that the new 
distribution is a non-central whose non-central value A is the power spectrum 

of the deterministic function describing the nonstationary signal. Finally, we test these 
results in the case of synthetic curves of gamma-ray bursts. We end up with a new 
formula for calculating the power spectrum uncertainties. This is crucial in the case 
of nonstationary short-lived processes affected by uncorrelated statistical noise, for 
which ensemble averaging does not make any physical sense. 
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1 INTRODUCTION 

The study of the temporal variability of time series in various branches of science and engineering has propelled the development 
of several techniques, both in frequency and time domains. Variability studies in the case of astronomical sources are crucial 
to gain insight over the dynamical and microphysical timescales, and therefore on the size of the emitting region as well as the 
nature itself of the emission process. This is of key importance in the X- and 7-ray domain, where remarkable flux variations 
are observed over timescales from days to ms. 

Fourier techniques are widely used in this field, as witnessed by the popular timing analysis package XRONofl 
IStella fc Angelinill 19921 ) included in the NASA heasoft packag^. The Fourier spectral analysis is fundamental in the study 
of stationary processes, since it provides an immediate physical interpretation as a power-frequency distribution. The Fourier 
power density spectrum (hereafter, PDS) in particular d ecomposes the total variance of a given time series to the different 
frequencies thanks to Parseval's theorem (e.g.. Ivan der Klis 1988: hereafter, K88). PDS analysis and related tools are suitable 
to both searching for possible periodic signals hidden in the data, and to characterising the so-called "red noise" connected 
with the presence of aperiodic variability. 

Practically, in the most general case one divides the time series in multiple adjacent intervals, over which the corresponding 
PDS is calculated. Finally, for each frequency bin the resulting PDS value and uncertainty are the mea n and standard deviation 
of the corresponding power distribution, as is routinely done by the dedicated XRONOS tool powspec (|Stella fc Angeliml992j '). 
However, the fundamental assumption behind this procedure is that each time interval represents a different sampling of the 
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same stochastic stationary process. The operat ion of replacin g ensemble averages with time averages of a single realisation 
makes sense only if the process is ergodic (e.g., |Priestlevl[l98ll ). 

Astronomical X-ray data rely on photon counting instruments. As such, measurement uncertainties are often dominated 
by the photon counting sta tistics, i.e. the Poi sson distribution. This translates into white noise in the PDS; adopting the 
normalisation introduced bv lLeahv et al.l (|l983l ) (hereafter L83), the power is known to follow a distribution with 2 degrees 
of freedom (xl). 

When additional, genuine variability of the source is also present, the resulting power distribution changes correspondingly. 
For instance, there is a class of X-ray sources, such as Cyg X-1, whose time series are the result of a process characterised by 
source-intrinsic correlated noise, for which time averages of PDSs from individual intervals are still meaning ful. These time 
series belong to a class of random processes whose PDS is still compatible with a rescaled xi distribution (K88: [Israel fc Stellal 
1 19961 ). 

In the case of highly nonstationary and short-lived events such as gamma-ray bursts (GRBs), the problem of a proper 
treatment of the Fourier PDS requires particular care. Dividing the light curve of a single GRB into several sub-intervals, 
and deriving an average PDS does not make any physical sense, given that the phenomenon is everything but stationary. 
Therefore, only a single sample PDS can be calculated over the entire observation duration. In these cases the statistical 
distribution of the PDS remains to be determined in the more general case of red deterministic variability due to the source. 
One has then to be careful with assuming the 100% uncertainties on the PDS provided by standard tools such as powspec: 
indeed, this relies on the X2 distribution, which does not hold any more in the more general case of a variable source with red 
noise. 

In the GRB literature, there have been different approaches. In one case, if one considers a set of different time series due 
to different GRBs as many realisations of the same stochastic process, then averaging the PDS of different GRBs still make 
sense. However, it must be pointed out that a strong assumption lies behind this: i.e., there is a unique stochastic process 
giving rise to the variety of observed GRB time profiles. This way one gai ns insight into the prope r ties of this general process, 
whose PDS is found to be described with a power-law, PDS oc f~^^^ l|Beloborodov et al.ll200ol : ISpada et al.ll200d ). In the 
other case, each GRB time profile is considered individually as the unique sample of a unique stochastic process, which is 
different from other GRBs. The statistics of its unique PDS is not known and is no more a X2- Monte Carlo simulations 
aimed at estimating the PDS u ncertainties by gener ating other (synthetic) samples of the same process cannot make use of 
the observed sample curve (e.g.. lUkwatta et al.ll201ll ). Indeed, such a procedure increases the noise variance and changes the 
xi nature itself of the power distribution. 

In this work we address the issue of a correct evaluation of the statistics of the PDS, and in particular of the power 
uncertainties, for a single sample of a nonstationary and short-lived signal, such as that of GRB time profiles. In this work 
time series are meant to be deterministic profiles affected by uncorrelated noise. We derived a formula for the variance of the 
PDS and show that it agrees with the known results of L83 in the pure white noise case. Finally, we test the validity of our 
results in the case of synthetic light curves of typical GRBs. 

Hereafter, time series are assumed to be discrete, equispaced, and with no data gaps. Our treatment assumes the noise 
to be purely statistical and uncorrelated, and we consider the two cases of Poisson or Gaussian statistics, suitable to photon 
counting detectors. PDSs are calculated assuming the Leahy normalisation (L83). We therefore neglect dead time effects, 
which are known to affect the statistics and suppress the variance of the resulting time series fe.g.. lMulleilll973l . K88). 



2 DESCRIPTION 

Let Xk (k — 0, . . . , N — 1) he a. time series observed within a time window with a duration T. The corresponding bin times are 
tk = kT/N . The observed series represents a sample function of the true time series we would have observed in the absence 
of statistical noise (e.g., with an infinite collecting area detector). Hereafter, random variables are written in bold. Each Xk 
is therefore a single sample of the random variable Xk, whose expected value and variance are defined as E{xk} ~ rjk, and 
E{{xk — Vk)^} ~ (^k- The random variables are assumed to be independent, E{xk xi} = rik rji {k 7^ I). 
The discrete Fourier transform (DFT) amplitudes are defined by 

a, = Y^xke--^-'- + (1) 

fe=0 

jv/2-1 

= ^ E fc = 0,...,iV-l, (2) 

j = -JV/2 

where the corresponding j-th frequency is /j = j /T. The highest frequency is the Nyquist frequency, /jv/2 
Leahy normalisation the power spectrum is defined by 

k,l 

or, equivalently. 



= 7V/2 r. In the 
(3) 
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= ^ J2 (^'^(^ " ').?/^) • (4) 

^'^ k.l 

= "^k-o "^ft' expected total variance. As shown in Section [Bl this happens to coincide with the total number of 

counts in the specific case of Poisson statistics, which explains the reason for the choice of this name. From equation ((3]) the 
expected value of the random variable Pj is derived straightaway, and is 

E{P,} = 2 + -^ ^^.r,,e-'('=-"^'"^ , (5) 

where we used E{xj.} — a^+rjl. The first term in the right-hand side of equation ((5]) accounts for the statistical noise variance, 
also called "white" because it does not depend on frequency. When the signal is constant, rik = rj (Vfc), Pj {j = 1, . . . , N/2 — 1) 
is known to be xi distributed, except for the Nyquist frequency, for which Pjv/2/2 is Xi distributed and must be treated 
separately (L83). In the pure white noise case the second term vanishes, thus leaving E{Pj} — 2. 

In general, the second term can also be seen as the DFT of the autocorrelation function (ACF) defined by 

JV-l 

Rj = ^'?fe')fc+. , (j = 0,...,7V-l), (6) 



assuming the periodic boundary condition rjk+N = ??fe. This is the discrete versio n of the Wiener-Khinch in theorem stating 
that the power spectrum of a time series is the Fourier transform of its ACF (e.g.. iPapoulis fc Pillaill2002t ). 
We define the power density spectrum of the deterministic function Pj^^ by 



AT-l 



2 

(7) 



' iVph I ^ 

so that equation (O can be written as 

E{P,} = 2 + 7^('" (8) 

So far no assumption was made on the kind of the distribution of the random variables Xk- Hereafter, we distinguish 
between the Gaussian and the Poisson noise cases. 



2.1 Gaussian noise 

Each random variable Xk is distributed according to a normal N{r]k,(Jk)- We can express equation Q as the result of matrix 
products: 

Pj = AX , Aki = ^ cos (27r(fc - l)j/N^ (9) 

where X is the column vector whose fc-th row element is Xk, and Aki is the ( fc, I) element of A. The matrix A is a positive- 
definite quadratic form. As such, from the algebra of the quadratic forms (e.g.. lEnnis fc Johnson|[l993l ) the distribution of the 
random variable Pj is a non-central xi{\) if and only if the two following conditions are fulfilled: 

A = ASA 

rank(A) = r ' ^ ' 

where r is the degrees of freedom, and 



, (11) 



with H being the column vector whose fc-th element is rjk- Equation follows from the definition of Pj''' given in equa- 
tion ((Tjl. E is the covariance matrix, so its (fc,/) element is given by 

Efei = CoY{xk,xi) = E{xkXi} - E{xk} E{xi} = Suicrl, (12) 

where 5ki is Kronecker's delta. To verify the first equation (jlOjl . we calculate the (fc, I) element of the matrix ASA, which is 
found to be: 

JV-l 

(A E A)ki = Aki + ^Y."''^ + ^ - 2m)j7iv) (^^ = 1, . . . , ^ _ (13) 



m=0 



The second term of the right-hand side is identically zero when all the variances cr^ are equal. More generally, the second term 
of the right-hand side of equation (|13|l is 0(1/A'ph) times Aki- Therefore, equation \1Q\ is approximately fulfilled in practical 
cases of interest. 

The rank r of A does not depend on H , so we can calculate it in the special case when rjk — rj (Vr;), i.e. the case of 
constant signal. We already know that Pj is xi distributed (L83), so it must be r = 2. 
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This proves that in the more general case of a varying signal described by a deterministic function whose discrete values 
are given hy r]k {k = 0, . . . , N and affected by pure uncorrelated Gaussian noise at, the Leahy-normalised power spectrum 
Pj distributes according to a non-central X2{^)j where A is P^^^ f equation 1111) . The expected value of Pj is already known 
from equation ((S)), while its variance ii^ 

Var(P,) = 2 (2 + 2 A) = 4 (l + ^ ^ r,, e^^'f^-"^/^) = 4 (l + P,*"') = 1, . . . , ^ _ l) (14) 



Vph 

k.l 



Appendix 1X1 reports the direct calculation of Var(Pj) and equation (|A7|) reports the exact formula for the variance. In the 
j = A''/2 case Pjv/2/2 satisfies the conditions (|10|) with r — 1, 

Var(P^/2) = 4Var(P^/2/2) = 8 (1 + 2 A) = 8 (l + -|- ^ r;, ,7, e^'C^-") = 8 (l -f P^%) , (15) 



Vph 

k,l 

where we used A = Pj^'/2/2. Equations p4ll5|) can also be written as a function of the expected power: 

Clearly, in the general case of uncorrelated noise and a significant power above the white noise level, i.e. when E{Pj} > 2, 
assuming a 100% uncertainty on the resulting power spectrum Pj, as assumed by widespread tools in X-ray astronomy such 
as XRONOS powspec overestimates the uncertainty. While the 100% uncertainty inherited from the xi distribution is correct 
under the assumption that all power is due to Poisson noise and int rinsic correl ated noise, in the case of uncorrelated noise 
here considered from equation (fTO)) the uncertainty is given by 2 E{Pj} — 1, and not E{Pj}. The two coincide only in 
the case of pure noise associated with a constant signal, so that E{Pj} = 2. The 100% uncertainty assumed by XRONOS 
remains correct in the presence of genuine stochastic variability due to correlated noise. 

In practice, when only a single sample series Xk is available of a short-lived process, such as the time profile of a GRB, 
the deterministic series r]k is not known a priori. However, one could constrain the probability density function (pdf) of Pj, 
i.e. the unknown A (equation [11} , with the only sampled value Pj. In principle, this makes no difference to taking the light 
curve with observed Xk counts as Xk ± y/xk instead of the unknown rjk ± ^/rjk for a Poisson process (the Gaussian case is 
formally the same). We estimate A, the best value for A, adopting a Bayesian approach: 

Vir,m) - ^-(^^^ , (17) 

where p{r,\\Pj) is the posterior function of the parameters (r. A) (r is fixed to either 1 or 2 depending on whether it is 
or not j — N/2) given the observed Pj. The likelihood function p{Pj\r,X) is merely the pdf of Pj, i.e. x^{\,Pj). The prior 
p(r. A) may include the knowledge one might have on A prior to measuring Pj, e.g. when a specific shape of the deterministic 
PDS is expected. In the most general case, we assume a uniform prior. The term p{Pj) normalises the likelihood function. 
For a given Pj A is chosen so as to maximise the posterior probability. In this case our approach is equivalent to a maximum 
likelihood estimation (MLE). In Appendix[C]we show that it is A = (Pj < 2), and that A rapidly converges to Pj {j < N/2) 
and to P]v/2/2 (j — N/2) for Pj > 2. We conservatively assumed A — Pj (= Pjv/2/2 for j = A'^/2) for all values of Pj, to 
avoid the risk of underestimating the variance. Therefore, for a single time series Xk consisting of an unknown deterministic 
function affected by uncorrelated white noise equations (|14p and (|15p are approximated to 



Interestingly, the case of a constant va,r iance (all are equal), for which the first equation (|10p is fulfilled exactly (see 
equation [13]), was discussed bv lGrothI (|l975h . Apart from a scale factor of 2 in the definition of power, the probability density 
function he derived is precisely that of a non-central chi square with 2n degrees of freedom, and a non-central parameter 
given by the deterministic (or "signal', as he called it) power (see equations 12 and 14 therein). His treatment considered the 
case where the power is the sum of n terms due to as many frequency bins, while here we consider the n = 1 case. 



2.2 Poisson noise 

Each random variable Xk is distributed according to a Poisson distribution with expected value r]k. Since Xk are no more 
normal, we cannot exploit the properties expressed by the conditions (|10p for the Gaussian case. In Appendix |B] we calculate 
the corresponding variance of Pj. Equation (|B4|) reports the exact formula for Var(Pj). 
The main results are the following: 

• if the observed counts Xk are so high as to ensure the Gaussian regime [xk 3> 1) the results are the same as those discussed 
in Section [2. II since the overall process essentially becomes Gaussian. 



3 If y is xlW distributed, then E{y} = r + A, and Var(y) = 2 (r -|- 2 A). 
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Figure 1. Top left: example of a synthetic curve of a FRED-shaped GRB light curve. The thick solid line represent the deterministic 
process. Top right: power spectrum of the sample curve shown in the top left panel. Uncertainties have been calculated following 
equation III8I I. The solid line shows the power spectrum of the deterministic function, while the dashed line shows the white noise level, 
2. Mid right: ratio between the uncertainties provided by the tool powspec and the correct values determined from the MC simulations 
for the corresponding frequency bin. Bottom right: ratio between the calculated uncertainties and the correct value from MC simulations. 
Bottom left: the noise-subtracted power of the same sample function has been binned up so as to ensure 3cr significance. The upper limit 
is at 3(7 . 



• If A'ph S> 1, even if the individual Xk are in the low-count regime, equations (|14H18p still hold, and in particular Pj still 
distributes according to a non-central xiW, ^ in Secti0n r2.ll 

• If A'^ph ~ few, equations (|14H18|) do not hold any more and the distribution of Pj deviates from a xiC-^)- 



3 APPLICATIONS: GAMMA-RAY BURST LIGHT CURVES 



We produced a number of GRB synthetic light curves to test the validity limits of the results we derived for nature of the 
the power spectrum distribution and summarised by equations H14H18[) . GRBs are particularly suitable to this aim, given 
their nature of highly nonstationar y and short-lived phenomena. We adopted the fast-rise exponential decay (FRED) profile 
as modelled bv lNorris et al.l (| 19961 ). This function satisfactorily describes the temporal behaviour of the simplest example of 
GRB light curve, which consists of a single pulse modelled as 



A exp 
A exp 



Ti ) 



(19) 



where tmax is the peak time, Tr and Td are the rise and decay times, respectively, A is the normalisation and p is the peakedness 
(when V = 1 the profile is a s imple exponential, when = 2 it is a Gaussian). For a typical FRED it is Tr/r^ < 1, with 
an average value of ~ 0.3-0.5 H orris et al] 1 19961 ). The continuous power density spectrum of a FRED with p — 1 (double 
exponential) can be calculated analytically, and apart from a normalisation term is found to be 



A^ {Tr + Td) 



(20) 



[l + {2TVjyTry] [l + {2 

Several examples of power spectra obtained for the more general case of a FRED with p 7^ 1 are discussed bv lLazzatil (|2002h 



We started with a FRED with the following parameters: Tr = 10 s, Td = 30 s, p = 1.5, A — 1000 counts bin 



= s. 



superposed to a constant detector background with an average intensity of 1000 counts bin~^. We generated 5 x 10'^ samples 
of this pulse with a bin time of 64 ms, assuming Poisson statistics. The total number of bins amounts to 4096 (= 2^'^). Given 
the large number of counts per bin, this is equivalent to the Gaussian case. Figure [1] displays the synthetic curve of the 
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Figure 2. Left: Leahy power distribution at u = 0.061 Hz derived from 5000 synthetic sample curves of the FRED pulse of Fig.[T] The 
expected distribution is shown with a solid line, and corresponds to a non-central x^C-**); with A = 36.41. Right: the same at i/ = 0.076 Hz, 
where A = 10.69. 



deterministic model as well as one out of the simulated samples. The deterministic PDS is used for comparison to assess the 
goodness of the PDS of the sample curve, in particular of the uncertainties derived for the power adopting equation (|18p as a 
function of frequency. This is shown by the ratio between the calculated a and the scatter of the power of the synthetic PDSs 
observed for each corresponding frequency bin (bottom right panel of Fig. [!}. Apparently, the ratio ranges between 0.5 and 
2. An overall comparison with the analogous ratio between the uncertainty provided by powspec and the corresponding value 
determined from the MC simulations (mid right panel of Fig[l| shows that the improvement is noteworthy. Furthermore, the 
accuracy of equation (|18p is evident when the PDS is dominated by the deterministic power of the signal (at low frequencies). 
At high frequencies, where the white statistical noise dominates, the power uncertainty is systematically overestimated up to 
a factor of 2. However, compared to the values provided by powspec often underestimated by up to an order of magnitude, 
it still represents a significant improvement, particularly in the more conservative direction of overestimating rather than 
underestimating uncertainties. 

So far this proves that equation ((18} overall provides a satisfactory means for estimating the uncertainty on the PDS 
of a sample curve. We go further and test whether the distribution of the individual Pj is indeed a non-central chi-square 
distribution. To this aim, from the PDS of the sample curve considered above we chose two frequencies and derived the 
corresponding power distribution from the synthetic PDSs. The result is displayed in Fig. [2] The expected non-central 
parameter for the corresponding X2{^) was calculated with equations (f71ll4p. 

Figure [3] shows the same FRED pulse 200 times fainter superposed to a correspondingly fainter constant background 
of 5 counts bin~^. The single Xk variables cannot be approximately assumed to be normally distributed, so that we may 
test the Poisson regime. Still, the total number of counts is still very large, iVph = 2.3 x 10*. This means that the results 
obtained for the Gaussian case should still hold. Indeed, both the comparison between the calculated uncertainties and the 
scatter of the corresponding power from the simulated PDSs gives similar results to the previous case (bottom right panel of 
Fig. [3]). Likewise, the power distributions of individual frequency bins are fully compatible with the corresponding expected 
non-central chi-squares, as shown by Fig. 



4 SUMMARY AND CONCLUSIONS 

We investigated the nature of the statistical distribution of the power density spectrum of a single, nonstationary, and short- 
lived time profile on a theoretical ground. This treatment assumes the time series to be deterministic profiles affected by 
uncorrelated noise. In other words, the time series here considered consist of a set of statistically independent, Poisson and 
normally distributed random variables, whose expected values represent the deterministic function of the varying signal to be 
studied. 

We demonstrated that the probability density function of the power is a non-central X2(-^)i whose non-central parameter 
A corresponds to the power of the deterministic function. This holds in the Gaussian case, as well as in the Poisson case, 
provided that the Gaussian limit of A'^ph ^ 1 is fulfilled (A^'ph being the total number of counts) . As a consequence, we provided 
a new formula for calculating the correct uncertainty of the power at each frequency as a function of the observed power itself. 
We finally showed the agreement with simulated light curves of typical GRB time profiles. These results provide a statistically 
solid basis to a proper treatment of power density spectra in the case of nonstationary and short-lived time series affected by 
uncorrelated noise. 
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Figure 3. Same as in Fig. [T] Here the FRED has the same profile, but it is 200 times less intense. The background level is proportionally 
lower, 5 counts bin~^. This example fits in the low-count rate Poisson regime (Afph = 2.3 X 10*). 
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Figure 4. Same as in Fig. [J] referred to the pulse of Fig. [S] 
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APPENDIX A: VARIANCE OF THE POWER (GAUSSIAN CASE) 

The central moments of a random variable Xk normally distributed as N{rik,o-k) are 

E{{xk - 7jkf'+'} = (Vi), E{{xk - Vkf} = '^l E{(xk - r^k)*} = 3at (Al) 

where the different variables at different values of k are meant to be independent, so E{xkXi} = rjkrii, {k ^ I). They can be 
used to calculate the corresponding noncentral moments through the following. 
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The noncentral moments are given by 

E{xk} = rik, E{xl} ^ i]l + al, E{xl} ^ i]l + Srjkal, E{x'i} = r]t + filial + 3at. (A3) 
The variance of Pj is calculated directly from equation (|3]): 



Var(P.) = (^) (£{P?}-i?'{i'.}) (A4) 



, 2 

2iTi{k-l)j/N 1 



k.l,m.n 



where we used the definition of A'^ph. The terms of equation (|A4[) must conveniently be separated based on the different kind 
of moments. To simplify the notation, we define cj = 2Tr j/N. 



Var(P,) = + 2 E E{xl}E{x,}(^e-'^''-'^ + e^'^'^-") + ^ (2 + e^-'*^-") + (A5) 

k ki^l kj^l 

k,l^m,n^^ k,l k,l^m,n 

r, 2 uiUk-l) 

k.l,m 

We adopted the following notation: ^ ^ ^ is a sum over k, I, and m, and is meant to exclude all the equality cases (fc 7^ /, 
k ^ m, I ^ m). Replacing the values of the corresponding moments, 

= E(6%-^ + 3-^) + 6 E {e-''-'' + e-"'<'=-") + Y^^^laf + vUf + nUl) (2 + e^'f^-") + (A6) 

k k^l kj^l 

E2 /, uii[l-m) I uii[2k-l-m) , -uii(2k — l — m)\ \ ^ 2 2 o \ ^ 2 uii(k — l) 

crkViVrn{4:e '+e '+e ) ~ / .""fc""' ~ ^ / . '^mVkVie 

k.l,m.^ k,l k,l.m 



After a few passages, by adding and subtracting the terms excluded in the sums in equation (|A6|) . one ends up with the 
following: 



Var(P,) = 4 (l + ^E'?^'?'^'^'*'""""')+ 

2 2 4ni(k-l)j/N , ST^ 2 / 2wi(2k-l-m)j /N , -2ni(2k-l-m)j / N \ 



" k.l k.l.m 



Equation (|A7|) is exact. For the Nyquist frequency, equation (|A7|) becomes 

Var(P^/2) = 8 (l + ^ E^*'?'^""*'"'' ) (^S) 



in agreement with equation p5|) . For j 7^ A/2 the second term in the right-hand side of equation (|A7|) is identically zero when 
all (Tfc = (T (Vfc), and it can be neglected when the (t^'s are comparable with each other, being 0(1/Aph) times the first term. 
Equation (|A7|) can be approximated by 



^ph 

in agreement with equations (|14I15|) 



APPENDIX B: VARIANCE OF THE POWER (POISSON CASE) 

The case of a process Xk affected by Poisson noise is a more general case than that of a Gaussian noise, since the latter 
corresponds to the former in the high count rate regime, i.e. when it is E{xk} = rjk ^ 1- For a Poisson process, the central 
moments are 
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E{{xk-rjkf}^rik, E{{xk - rikf} ^ Vk, E{{xk - ritf} = T^k + Sr,l (Bl) 
From equation (|A2|) the corresponding noncentral moments are 

E{xk}^m, E{xl} = rjl+7^k, E{xl} ^7^l + 34+m, E{xt} ^ r,t + 6r,l + 7ril + 3^- (B2) 

In the Poisson case the normahsation constant A'^ph — X/^o^ ^® expected total counts. The definition of the power 
spectrum Pj is the same as that of equation (|3]). Using the first two moments of equation (|B2|) . calculating the expected value 
of Pj is straightforward, and is found to be the same as equation ((5|. 

Calculating the variance of Pj in the Poisson case is formally the same as the Gaussian case up to equation (|A5|) . i.e., 
prior to substituting the specific values of the moments. At this point the two cases must be treated separately. Replacing the 
moments of equation (|B2|) in (IA5|) . and defining co = 2n j/N as before, it becomes 



(^) Var(P,) = Y.^6vl + 7d + Vk) + Y.(6vl+'2Vk)vi{e-^^'~'^+e--'^''-'^'^+Y^ (B3) 

k k^l k^l 

/ I 2 I 2\ 2ui(k — l) I ^ / A uji(l—m) , ijji(2k — l — m) , —uji(2k — l — m) \ 

{■nkVi + VkVk + VkVi )e \+ VkViVm[4:e '+e ^ '+e 'j 



k,l,m,j^ 



uji(l — m) 

■ilkViVme 



Similarly to what was done in section [A] adding and subtracting the excluded terms in the sums, after a few passages one 
ends up with 



k,l 

4 [V^ 4,wi(k-l)j/N . I 2wi(2k-l-m)j/N , -2wi(2k-l-m)j / N 

j:j2- [Z^VkVie +2_^VkViVm{e +e 



As done for the Gaussian case, we have to treat the Nyquist frequency separately. When j — N/2, equation (|B4|) becomes 
Var(P./.) ^ 4 (2 + J-) + A. (2 + 2 ^ ^ ^ ^^^^ 



which is equivalent to equation HA8|) in the Aph ^ 1 limit. In the special case of a constant signal rit = rj (Vfc), equation (|B4|I 
reduces to 

in agreement with the results of L83. In the more general case of a nonstationary signal rjk, in the limit Aph ^ 1, equation (IB4I) 
can be approximated by 

Var(P,) ^ 4{l + -l-J2vkVie'-'''-'^'^'')+ (B7) 
^'^ fc,i 



k.l k,l.m 



Not surprisingly, equation (|B7[) is the same as (|A7[) upon replacing with rik, as expected for a Poisson variable in the 
Gaussian limit. As discussed in Appendix [X] the result in equation (|B7]) can be approximated by equation (IA9|) . provided 
that Aph 3> 1. When the Gaussian limit is not satisfied, i.e., when Aph is just a few, we note that the relation between 
expected value and variance for a noncentral chi-square distributed random variable with r — 2 degrees of freedom (j 7^ A/2) 
is not fulfilled: 

E{P,} = 2 + ^ yvkVie^^'^"'''"^'' = r + X (B8) 

JVph ^ — ' 
k.l 

Var(P,) = 2{2+J-Y Vkme'^'''-'^^^'') + ^ (l + ]^ E VkVie'^'"-'''''') + (B9) 

k.l k.l 



k,l k.l,m 

= (r + 2A) + -^(r + 4A) + -^[...] / 2 (r + 2A). 
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We conclude that for A'ph ~ few, the distribution of the power spectrum at a given frequency is not a noncentral xi{\), 
as found in the Gaussian limit, and the expression for the variance to be used is given by equation (|B9|) . In the same regime 
of A^ph ~ few, Pjv/2/2 also deviates from a noncentral Xi('^) distribution, given that equation (jB5|) does not fulfil any more 
the corresponding relation between expected value and variance. The exact formula for the variance of the Nyquist frequency 
power therefore remains equation ([BSp . 



APPENDIX C: MAXIMUM LIKELIHOOD ESTIMATION 



The j < N/2 and j — N/2 cases must be treated separately, given that the probability density functions of the corresponding 
random variables Pj are non-central chi squares with r = 2 and r = 1 degrees of freedom, respectively. The purpose is to find 
the value for the non-central parameter A which maximises the likelihood function Xr(-^i Pj) for given measured value Pj. 



First, let us consider the j < N/2 case, for which it is r = 2. It is 



1 n 



^ /A cos 6 7/1 

eV J do 



(CI) 



When Pj = equation (|C1[I reduces to a simple exponential, so that A = 0. For Pj ^ 2, it is A = 0. For Pj > 2, A is found by 
requiring 



dX 



equivalent to 



= 



^ ^ hiVXPj) 

Pj U^Tp,) 



(C2) 



(C3) 



where 7„ is the modified Bessel function of the first kind and index n. At Pj ^ 1 the solution is A = Pj . Numerical solutions 
to equation IC3I are displayed in Figure [CTI (solid line), which shows how rapidly A converges to Pj as a function of Pj. 

In the j = A'^/2 case it is r = 1 and the interested random variable is Pjv/2/2. To simplify the notation, let us define 

X = Pjv/2/2, so it is 

1 



Xi(A,x) = 



^/2lfl 



cosh(%/^) 



(C4) 



When a; < 1 equation (IC4|) monotonically decreases for A > 0, so it is A = 0. When a; > 1, the maximum is found 
analogously to equation (|C2|I . thus 



dxl{\x) 



dX 







equivalent to 

1+ y/%/x 
1 - ^/fjx 



(C5) 



(C6) 



Analogously to the j < N/2 case, the solution is A < x and rapidly converges to A = x, as shown in Figure ICTl (dashed line). 

Summing up, in both cases it is A = for Pj ^ 2, and for Pj > 2 it asymptotically tends to Pj (P]v/2/2) for j < N/2 
(j = N/2). 

In the process of estimating the variance of Pj, we conservatively assume A = Pj {j < N/2), and A = Pjv/2/2 {j — N/2), 

so 



A : 



N/2 



/2 



(j < A^/2) 
U = N/2) 



(C7) 



By replacing equation (|C7|) into equations (fTll [T5|) equation (fTS)) is obtained. 
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